library(dplyr)
library(readr)
library(sf)
library(ggplot2)
library(mapview)Creating maps
Questions
- How do we create maps using R?
Objectives
- Learn how to plot iNaturalist observations on a map
- Learn how to create static maps with ggplot2
- Learn how to create interactive maps with mapview
Geographic concepts
Geographic data is data that has a location. There are various file formats for geographic data. Shape files for GIS applications, KML for Google maps, geojson for web applications.
Earth is a 3D sphere. Maps are 2D representation of a 3D sphere. Map projections are ways to represent a sphere as a flat surface. A coordinate reference system (CRS) defines the two-dimensional, projected map in your GIS relates to real places on the earth.
You can get boundaries for countries, states, cities, etc from various sources. I googled “Los Angeles county boundary shape” which had a link to “County Boundary | City of Los Angeles Hub - LA GeoHub” https://geohub.lacity.org/datasets/10f1e37c065347e693cf4e8ee753c09b I downloaded the shapefile for LA county.
You can also create your boundaries using GIS applications or GIS web applications.
Mapping iNaturalist data
iNaturalist data includes latitude and longitude, which means we can put the observations in a map. There are several packages to create maps. We will use ggplot and mapview packages.
loading R packages
Read data from the cleaned iNaturalist observation file.
inat <- read_csv('data/cleaned/cnc-los-angeles-observations.csv')Use names to see all the column names. “latitude” and “longitude” are the column names we need.
names(inat) [1] "id" "observed_on_string"
[3] "observed_on" "time_observed_at"
[5] "time_zone" "user_id"
[7] "user_login" "user_name"
[9] "created_at" "updated_at"
[11] "quality_grade" "license"
[13] "url" "image_url"
[15] "sound_url" "tag_list"
[17] "description" "num_identification_agreements"
[19] "num_identification_disagreements" "captive_cultivated"
[21] "oauth_application_id" "place_guess"
[23] "latitude" "longitude"
[25] "positional_accuracy" "private_place_guess"
[27] "private_latitude" "private_longitude"
[29] "public_positional_accuracy" "geoprivacy"
[31] "taxon_geoprivacy" "coordinates_obscured"
[33] "positioning_method" "positioning_device"
[35] "species_guess" "scientific_name"
[37] "common_name" "iconic_taxon_name"
[39] "taxon_id"
We use the sf package to add geographic data to our dataframe. st_as_sf() from sf package will take the longitude and latitude and add a geometry column that we can use for mapping.
- We pass in longitude and latitude columns to coors argument. Must wrap longitude and latitude in quotes.
- crs is coordinate reference system. 4326 is code for the WGS 84 CRS. It used by many GIS applications.
remove=FALSEwill keep the longitude and latitude columns in the dataframe
temp <- inat %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 4326, remove=FALSE)use names() to get a list of all the columns. A geometry column was added.
names(temp) [1] "id" "observed_on_string"
[3] "observed_on" "time_observed_at"
[5] "time_zone" "user_id"
[7] "user_login" "user_name"
[9] "created_at" "updated_at"
[11] "quality_grade" "license"
[13] "url" "image_url"
[15] "sound_url" "tag_list"
[17] "description" "num_identification_agreements"
[19] "num_identification_disagreements" "captive_cultivated"
[21] "oauth_application_id" "place_guess"
[23] "latitude" "longitude"
[25] "positional_accuracy" "private_place_guess"
[27] "private_latitude" "private_longitude"
[29] "public_positional_accuracy" "geoprivacy"
[31] "taxon_geoprivacy" "coordinates_obscured"
[33] "positioning_method" "positioning_device"
[35] "species_guess" "scientific_name"
[37] "common_name" "iconic_taxon_name"
[39] "taxon_id" "geometry"
use select to pick which columns to use.
inat_map <- inat %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 4326, remove=FALSE) %>%
select(id, user_login, common_name, scientific_name, observed_on, url, longitude, latitude, geometry) Use dim() to show the number of rows and columns. There are over 169K rows.
dim(inat_map)[1] 169045 9
static map
Use ggplot to create static map for the observations. geom_sf will use geometry column to produce a map.
ggplot() +
geom_sf(data = inat_map) 
interactive map
use mapview package to create interactive maps where you can zoom in and out.
Since there are over 169K rows, an interactive map will be very slow. I suggest not using mapview if there are lots of rows.
To speed up the interactive map, let’s filter the list of observations. Get all observations for Western Fence Lizard.
inat_lizard <- inat_map %>%
filter(common_name == 'Western Fence Lizard')Use dim to get number of rows. About 3000 rows.
dim(inat_lizard)[1] 2956 9
Create interactive map. You can zoom in and out. Click on observation to see the info.
mapview(inat_lizard)working with other geographic files
Let’s add LA county boundaries to the map.
I downloaded the LA county boundaries from https://geohub.lacity.org/datasets/lacounty::county-boundaries/explore
use read_sf() from sf package to read the shape file.
la_county <- read_sf('data/raw/County_Boundary/County_Boundary.shp')add LA County to maps.
ggplot() +
geom_sf(data = la_county) +
geom_sf(data = inat_lizard) 
mapview(la_county) +
mapview(inat_lizard) Exercise 1
Create a map for one species. Include the boundaries for LA County.
inat_finch <- inat_map %>%
filter(common_name == 'House Finch')
mapview(la_county) +
mapview(inat_finch)Exploring iNaturlist data in specific area
Lets look for all iNaturalist observations made in Exposition Park.
I used this Draw map boundaries to draw and download the boundaries of Exposition Park.
expo_park <- st_read('data/raw/boundaries_expo_park_area.geojson') %>%
st_transform(4326)Reading layer `boundaries_expo_park_area' from data source
`/Users/wyk/Development/science/city_nature_challenge/NHMLA_workshop/CNC-coding-workshop_quarto/lessons/data/raw/boundaries_expo_park_area.geojson'
using driver `GeoJSON'
Simple feature collection with 1 feature and 2 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -118.2915 ymin: 34.01096 xmax: -118.2829 ymax: 34.01806
Geodetic CRS: WGS 84
plot map of Expo Park.
ggplot() +
geom_sf(data = expo_park) 
mapview(expo_park) We want to get observation inside Expo Park.
You should check if the crs for the inaturalist data and the Expo Park are the same
st_crs(expo_park) == st_crs(inat_map)[1] TRUE
Use st_intersection() to get all observations that inside of Exposition Park.
inat_expo <- inat_map %>% st_intersection(expo_park)Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Use dim to get row and column count. 169K in LA county. 2500 observation in Expo Park.
dim(inat_map)[1] 169045 9
dim(inat_expo)[1] 2557 11
Create map of all observations in Expo Park.
ggplot() +
geom_sf(data = expo_park) +
geom_sf(data = inat_expo) 
mapview(expo_park) +
mapview(inat_expo) Use alpha.regions to set opacity. use col.regions to set color.
mapview(expo_park, alpha.regions=0.3, col.regions="#333333") +
mapview(inat_expo) Exercise 2
Create a map for one species that are inside of a specific area
- Used Draw map boundaries to draw and download an area that you are interested in.
- Save the file to the
data/rawdirectory. - use
st_readto read your boundary data. - use
st_intersectionto get observations inside a boundary :::::: answer
expo_park <- st_read('data/raw/boundaries_expo_park_area.geojson') %>%
st_transform(4326)Reading layer `boundaries_expo_park_area' from data source
`/Users/wyk/Development/science/city_nature_challenge/NHMLA_workshop/CNC-coding-workshop_quarto/lessons/data/raw/boundaries_expo_park_area.geojson'
using driver `GeoJSON'
Simple feature collection with 1 feature and 2 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -118.2915 ymin: 34.01096 xmax: -118.2829 ymax: 34.01806
Geodetic CRS: WGS 84
inat_finch <- inat_map %>%
filter(common_name == 'House Finch') %>%
st_intersection(expo_park)Warning: attribute variables are assumed to be spatially constant throughout
all geometries
mapview(expo_park) +
mapview(inat_finch):::
Exporting maps
We can export the static maps created with ggplot. We can not export the interactive maps created with mapview.
Assign the map to an object. Then run ggsave() to save our map The first argument we give is the path to the file we want to save, including the correct file extension. You can save as jpb, pdf, tiff, png. Next, we tell it the name of the plot object we want to save. We can also specify things like the width and height of the plot in inches.
my_map <- ggplot() +
geom_sf(data = expo_park) +
geom_sf(data = inat_expo)
my_map
ggsave(filename = 'results/expo_park_observations.jpg', plot = my_map, height = 6, width = 8)